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Abstract 

Most of the physically based techniques for rendering translucent ob- 
jects use the diffusion theory of light scattering in turbid media. The 
widely used dipole diffusion model JMLH01 applies the diffusion-theory 
formula derived for the planar surface to objects of arbitrary shapes. This 
paper presents first results of our investigation of how surface curvature 
affects the diffuse reflectance from translucent materials. 

1 Introduction 

Translucent materials, such as human skin, marble, wax, fruits, more scatter 
light than absorb it. Therefore, when a photon enters such a material, it un- 
dergoes many scattering events under the surface before it leaves the material. 
Such a light behavior is well described by the Bidirectional Surface Scattering 
Distribution Function (BSSRDF) [NRH+77] . Based on the light diffusion the- 
ory, Jensen et al. [JMLH01] suggested the dipole diffusion model for BSSRDF. 
This model applies an expression for reflectance from a turbid half-space to 
arbitrarily shaped objects. The multipole [DJ051 IDJ06) and quadpole [DJ08J 
models have been suggested to describe more complicated geometries - a mul- 
tilayered slab (or half-space) and a right-angle corner, respectively. Jensen et 
al. |DJ08j showed that a big variety of shapes can be rendered by combining 
photon tracing and a scheme for interpolating between dipole and quadpole and 
between quadpole and multipole models wherever appropriate. However, they 
do not focus on how the BSSDRF itself changes as a flat surface is replaced with 
a curved one. It is difficult to devise how their interpolation scheme can be used 
with approaches that do not use photon tracing - for example, the curvature- 
based method [Kol07 . Our goal is to investigate how inclusion of curvature may 
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change the diffusion BSSRDF model. A BSSRDF model that includes curvature 
effects could be easily incorporated into many existing approaches for rendering 
translucent materials. We present here preliminary results of our study. 



2 Diffusion Equation 

Under the assumption that light scattering in a turbid medium dominates ab- 
sorption, light transport in it is well described with the diffusion theory [Far92] . 
The fluence rate 4 , (r) obeys the modified Helmholtz equation [F ar92| 

A* - of r * = -D- X *(r - r ) (1) 

where <7tr — \Z3&a{&'s + a a) is the effective transport coefficient, a' s is the re- 
duced scattering coefficient, a a is the absorption coefficient, D — ^ a , 1 +c , s is the 
diffusion coefficient. We refer the reader to [JMLH011 [Far92j for explanation of 
the physical meaning of the quantities. In the above equation, we assume that 
there is a single source in the medium, and it is located at a point ro. 

Let us first consider the case of translucent material occupying the half-space 
z > 0. The point source is at ro = (0, 0, zq). Farrell et al. [Far92] showed that 
quite an accurate solution can be obtained by using the boundary condition 
^\z=-z b = and putting the image source at the point ro = (0,0, —zq — 2zf>), 
where zt, = 2 AD, and A is calculated as described in [JMLHOll [Far92j. The 
resulting fluence is 
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where n and r^. are the distances to the source and image source, respectively; 
that is, 



n = [ {z - Zo f + p ^ (2) 

r 2 = [(z + z + 2z fc ) 2 +p 2 ] 1/2 (3) 
The reflectance is calculated from the fluence using the formula 

R = -L>W (4) 
where the gradient is evaluated at the interface. In the planar case, this gives 



R(p,z Q ) = — 
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+ {z + 2z b )(a tr + -t) 
where n and r 2 are calculated for z = 0. 



(5) 



The dipole diffusion model [JMLHOlj applies the above formula to an arbi- 
trary shaped air-material interface by calculating r\ and r 2 as the distance from 
a point being shaded to the source and image source, respectively. 



3 Exact Solution for a Sphere 

Suppose the turbid medium is confined within a sphere having the radius Rq and 
the center at z = Ro- In addition to Cartesian coordinates, we will also use the 
polar system of coordinates with r counted from the sphere center and 6>counted 
from the z axis. We assume that i?ois much bigger than the mean free path for 
photons scattered in the medium, we can use the same boundary condition as in 
the planar case - namely, the fluence rate vanishes at a distance of Zb from the 
sphere surface. In other words, <I> is zero at a sphere of the radius R = Ro + z&. 
We will solve eq. © with the boundary condition x P| r= i?, = following the 
method described in |Mat71| . The solution of the modified Helmholtz equation 
Ovith the zero boundary condition on the sphere r — R can be written as 
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where r' is the distance of the point source from the sphere center; that is, we 
suppose thatro has the polar coordinates r = r' and 9 — 0. The functions I v (r) 
and K u (r) are the modified Bessel functions [M A70] . The constants A m and B m 
are determined by stitching the solutions [B] at the sphere r = r' . The function 
$ is continuous, but its derivative is not. In a manner similar to that used in 
|Mat71j . we integrate eq. [9] over an infinitisemally thin region confined by parts 
of spherical surfaces with radiuses r = r' + e and r = r' — e and containing the 
pointro. We utilize the Gauss theorem and get 









\~dr 


r'+e &r 


r' — e 7 



1 



(7) 



where is the solid angle variable. The delta function 6(d) can be decomposed 
in terms of the Legendre polynomials as [M at71j 
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Substituting eq. © into eq. ([7]) and calculating the derivatives from eq. 
([6]), we arrive at an equation for yl m and B m . One more equation for them is 
obtained by requiring continuity of ^ at r = r' . Solving the resulting system of 
two equations, we get 
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where 



r = (r 2 +r' 2 - 2rr' cos (9) 1/2 , 

and we used equality 10.2.35 from [MA70J. 

To find the reflectance, we choose r' = Rq — zq, apply eq. |4]and set r = Rq 
and get 
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4 Results 

We calculated the reflectance from translucent spheres of various radiuses. The 
incident light is a pencil beam entering a sphere at x = 0,y=0 . Ideally, we 
should consider a line of sources situated along the z axis. But it was shown in 
|Far92| that they all can be replaced with a single source located z = 1 /{a' a +a a ). 
The plot below shows how the reflectance depends on the distance from the point 
of light entrance measured along the surface (that is, the length of a geodesic 
connecting the entrance point and the point of interest). The calculations were 
done for the scattering coefficient a' s = fr7im _1 and absorption coefficient a a — 
O.OlTOm -1 (note that in |Far92) . the same quantities are designated as fi' s and 
fx a , respectively). These values of the scattering and absorption coefficients 
are typical for human tissue (see [JMLHOl] ). It can be seen that in this case, 
the difference between the exactly computed reflectance and that found by the 
dipole diffusion model becomes noticable only when the radius approaches 1 cm. 

Figure 1 above shows visualization of light reflection from spheres having a 
radius of 1 cm in two cases - a potato, on the left, and marble, on the right. As 
for the plot given below, we assume that a sphere is lit up by a stencil beam 
entering the sphere at the center of the image. The left part of each of the 
image corresponds to the exact calculation we describe above. The right part 
is computed using the diffuse dipole approximation. We used the measured 
values a' s and a a reported in [JMLHOf] , Because the amount of reflected light 




Figure 1: A spherical potato (left) and a marble sphere (right) illuminated with 
a stencil beam, which enters at the image center, normally to the image plane. 
Each of the spheres is rendered using the exact solution proposed (left part of 
a sphere) and the dipole diffusion model (right part of a sphere). 

decays with distance from the entrance point very rapidly, we applied the tone 
mapping operator to a calculated HDR image. We chose the logarithmic map- 
ping operator [DMAC03], as it is simple and robust, and a source code for its 
implementation is available on the web. 

Probably, we could anticipate in advance that the diffuse dipole model would 
underestimate the reflectance. However, our investigation shows that this under- 
estimation is small when curvature radiuses are of the scale of several centimeters 
and more for such materials as marble, potato, human tissue. 




Distance along surface (mm) 



5 Future Work 



The investigation presented here definitely lacks comparison of analytical re- 
sults with Monte- Carlo simulations. We are working on this and plan to report 
them elsewhere when the work is complete. Also, we would like to consider 
the case of arbatrarily curved surfaces. It would be interesting to try to build 
a phenomenological model for reflectance from a translucent material with an 
arbitrary surface. It can be sought as a function of principal curvatures at 
the point of light entrance. An approximate solution for slightly curved sur- 
faces (Appendix A presents an approximate solution for the modified Helmholtz 
equation in the case of a curved boundary) can serve as a base in attempts to 
construct a phenomenological model. Monte-Carlo simulations can be used for 
validation of such a model. A big potential of the phenomenological approach 
to constructing BSSRDF models has been proven by successfull development of 
an empirical BSSRDF model described in jDLR+09| . A BSSRDF model includ- 
ing surface curvature could be incorporated into the curvature-based method 
|Kol07| . It could be used for investigating perceptional effects, such as color 
shift at the terminator line |GreQ4] . 
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A An Approximate Solution for Slightly Curved 
Surfaces 

We consider the modified Helmholtz equation 

A* - cr t 2 ,,* = -D^Sir - r ) (9) 

where r is a point in three-dimensional space, and r = (x, y, z) in the Cartesian 
coordinate system. The source is located at the point r = (0,0, a). We are 
searching for a solution of this equation for z > f(x,y) with the boundary 
condition 



*!*=/(*,„) = o 



where z — f{x, y) is the equation describing some surface. Let's assume devia- 
tion of the surface from a plane to be small and introduce a small parameter, e, 
to describe this, so that the surface equation becomes 

z = ef(x,y), 

Let us search a solution to the modified Helmhotz equation as a series in e. 

* = *o + e*i +e 2 *2 + ••• 
Then *S? is the solution for the planar case, which is given by 

*o = 1 
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where r\ and r 2 are the distances to the source and image source, respectively; 
that is, 



n - [(z-a) 2 +x 2 +y 2 ] 1 / 2 
r 2 = [( z + a) 2 +x 2 +y 2 } 1 / 2 



(10) 
(11) 



We know the Green function for the planar case, so let's try to find the 
boundary condition for at the plane z — 0. It can be found by extending 
in z. We know that * takes a zero value at the surface; therefore, we have 



*U=e/(x,j/) = = *o|«=o 
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Each term in the e series should be equal to zero; therefore, we get 
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Therefore, 'I'ifor z > can be found as 
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where G is the Green function for the Helmholtz equation and boundary con- 
ditions of the form "J/| z= o = tp(x, y). Namely, 



G(x,y,z,x',y',z') 
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where r\ and r 2 are the distances to the source, (x, y, z), and the image source, 
(x,y, — z), respectively; that is, 



r[ = [(x' -x) 2 + (j/ -y f + {z' - z?] 1 '* (14) 
r' 2 = [(x' - xf + (y' - yf + (z' + zf] 1 ' 2 (15) 



We will calculate the integral in eq. [T3] using the Fourier transform. Repeating 
the reasoning of |Fei61| . we find that 
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and 
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Without loss of generality we may suppose that that the entrance point is 
(0,0,0), and the (x, y) plane is tangent to the object surface at the entrance 
point (this can always be achieved by a linear change of space coordinates). We 
assume that the interface is a smooth surface. It can then be locally represented 
as 

* = f(x,y), 

where f(x,y) is a smooth function, and its Taylor series expansion begins from 
quadratic terms. Choosing the x and y axes along the principal curvature 
directions, we get the following local surface representation 



f{x,y) = ^ (k lX 2 + k 2 y 2 ) 



(18) 



where ki,k 2 are the principal curvatures. The above equation is accurate to the 
second order in x, y. 

After substituting [TBI into [T3l we come to expressions kix' 2 q^- ' z ^ | ,_ 
and kny' 2 ^ ^^ ' z - | z /_ inside the integral. Their Fourier transforms can be 

calculated by taking the second derivative of the Fourier image of ^-^gjr ' z | z'=o 
in 1 171 with respect to q x and q yi rescpectively. Thus, the contribution to "Ji com- 
ing from the x coordinate is (the y part can be obtained by replacing x with 

y) 
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Integration by parts gives 
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Taking into account that 
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we arrive at 



where 
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Adding the contribution coming from the y coordinate, we get 
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